cd  "E:\Replication Data for Media Freedom and Public Trust in Science" 

use Data_Gu_Main, clear

**# 前期准备
***生成控制变量的宏
*个人层面： female  emp  income 
global cov ( i.female i.emp i.income )
*国家层面： loggdppc logpop  log_nobel_pc log_qs500_pc e_peaveduc e_pelifeex
global agg ( c.loggdppc c.logpop c.log_nobel_pc c.log_qs500_pc c.e_pelifeex )
**#所有控制变量
global control ( edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex )

***画图设置
set scheme s2color
grstyle init
grstyle color background white
grstyle set plain,horizontal grid
grstyle set color Set1
grstyle color ci_area gs12%40

**#Baseline Results
**#Table 1. Multilevel Model: Education Level
eststo clear

eststo h1: mixed y2 || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100 

eststo h2: mixed y2 edu female emp income || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100 

local var Mexp_conVDem_14after
eststo h4: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100

local var Mexp_binVDem_14after
eststo h5: mixed y2 edu female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex `var'_std || country:, mle
estat icc
estadd scalar ICC = r(icc2)*100

estadd scalar N2 = 108 , replace : h*

esttab h* using T1.rtf , b(3) se(3) nonote nobaselevel noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) nomtitles ///
refcat(female "{\i Individual-Level Variables}" loggdppc "{\i Country-Level Variables}" , nolabel) /// 
keep( _cons female emp income edu loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex Mexp_conVDem_14after_std Mexp_binVDem_14after_std ) ///
order(_cons)  ///
varlabel(_cons "Constant" female "Female" emp "Employment status" income "Income" edu "Education level" loggdppc "Log GDP per capita" logpop "Log population" log_nobel_pc "Log number of Nobel Prize winners per capita" log_qs500_pc "Log number of QS 500 universities per capita" e_pelifeex "Life expectancy" Mexp_conVDem_14after_std "Continuous media freedom exposure (V-Dem)" Mexp_binVDem_14after_std "Binary media freedom exposure (V-Dem)" ) ///
mgroup("{\i DV = Trust in Science (IRT)}") ///
title({\b Multilevel Model}) ///
stats( ICC N2 N, labels( "ICC (%)" "Country-level Observations" "Individual-level Observations") fmt( 2 0 0 ) ) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test).")

**#Figure 1. Multilevel Model: Marginal Effect of Media Freedom Exposure on Trust in Science as Moderated by Education Level
local var Mexp_conVDem_14after_std
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.edu || country:, mle
margins, dydx(`var') at(edu=(1 2 3))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("Continuous Media Freedom Exposure (V-Dem)") ytitle("Magrginal Effect of Media Freedom Exposure") xtitle("Moderator: Education Level") xsize(6) ysize(4) yscale(titlegap(-6))
graph save HLM_mod1.gph, replace 

local var Mexp_binVDem_14after_std
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.edu || country:, mle
margins, dydx(`var') at(edu=(1 2 3))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("Binary Media Freedom Exposure (V-Dem)") ytitle("Magrginal Effect of Media Freedom Exposure") xtitle("Moderator: Education Level") xsize(6) ysize(4) yscale(titlegap(-6))
graph save HLM_mod2.gph, replace 

*合并两个图
graph combine HLM_mod1.gph HLM_mod2.gph, ///
xsize(12) ysize(4) cols(2)
graph export "F1.png", replace

**#Table 2. Fixed Effects Model: Media Freedom Exposure and Education Level
eststo clear

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu     ,cluster(country)  keepsing
local r2=e(r2_a)
eststo m1: expint y2  `var'  beta
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu     ,a(country ) cluster(country)  keepsing
local r2=e(r2_a)
eststo m2: expint y2  `var'  beta
estadd local rfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_conVDem_14after
reghdfe y2     ///
c.`var'##i.edu $agg##i.edu ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m3: expint y2  `var'  beta
estadd local rfe "√"
estadd local bfe "√"
estadd local ctrl1 "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_conVDem_14after
qui reghdfe y2     ///
c.`var'##i.edu  $agg##i.edu $cov##c.`var'  ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m4: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

local var Mexp_binVDem_14after
reghdfe y2     ///
c.`var'##i.edu $agg##i.edu $cov##c.`var' ,a(country i.birthyear) cluster(country)  keepsing
local r2=e(r2_a)
eststo m5: expint y2  `var'  beta
estadd local ctrl2 "√"
estadd local ctrl1 "√"
estadd local rfe "√"
estadd local bfe "√"
estadd scalar r2a=`r2'
lincom primary-secondary
lincom primary-tertiary
lincom secondary-tertiary

esttab m* using T2.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) refcat(primary "{\i Effect of cumulative exposure to media freedom on}" , nolabel) /// 
varlabel(primary "Primary education or below" secondary "Secondary education" tertiary "College education or above") ///
mgroup("{\i DV = Trust in Science (IRT)}", pattern(1 0 0 0 0 0 0) ) ///
mtitles("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Binary Media Freedom Exposure (V-Dem)}") ///
stats(rfe bfe ctrl1 ctrl2 r2a N, labels("Country FE" "Birth year FE" "Individual-level controls × Media freedom" "Country-level controls × Education level" "Adjusted R²"  "Observations") fmt( 0 0 0 0 3 0 ) ) ///
title({\b Fixed Effects Model}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Further Analysis
**#Table 3. Mediation Effect of Political Trust
eststo clear

local var Mexp_conVDem_14after_std
local mv trust_govt
eststo m1: reghdfe y2     ///
`var'Xedu_nc     $cov##c.`var'    $agg##i.edu_nc  if `mv'!=.     ,a(country   i.birthyear) cluster(country)  keepsing
local total = _b[`var'Xedu_nc]
eststo m2: reghdfe `mv'     ///
`var'Xedu_nc  $cov##c.`var'   $agg##i.edu_nc  if `mv'!=.   ,a(country  i.birthyear) cluster(country)  keepsing
local im = _b[`var'Xedu_nc]  // stored for computing mediation effect
matrix SE = e(V)
local se_im = sqrt(SE[1,1])
eststo m3: reghdfe y2     ///
`var'Xedu_nc  `mv'  $cov##c.`var'   $agg##i.edu_nc    if `mv'!=.   ,a(country   i.birthyear) cluster(country) keepsing
local md = _b[`mv']
matrix SE = e(V)
local se_md = sqrt(SE[2,2])
estadd scalar med_p = `md'*`im'/(`total')*100
estadd scalar sobel_z = (`im' * `md') / sqrt((`im'^2 * `se_md'^2) + (`md'^2 * `se_im'^2))

local var Mexp_FOTPs_14after_std
local mv trust_govt
eststo m4: reghdfe y2     ///
`var'Xedu_nc     $cov##c.`var'    $agg##i.edu   if `mv'!=.   ,a(country   i.birthyear) cluster(country)  keepsing
local total = _b[`var'Xedu_nc]
eststo m5: reghdfe `mv'     ///
`var'Xedu_nc  $cov##c.`var' $agg##i.edu_nc     if `mv'!=.   ,a(country  i.birthyear) cluster(country)  keepsing
local im = _b[`var'Xedu_nc]  // stored for computing mediation effect
matrix SE = e(V)
local se_im = sqrt(SE[1,1])
eststo m6: reghdfe y2     ///
`var'Xedu_nc  `mv'  $cov##c.`var'   $agg##i.edu_nc    if `mv'!=.   ,a(country   i.birthyear) cluster(country) keepsing
local md = _b[`mv']
matrix SE = e(V)
local se_md = sqrt(SE[2,2])
estadd scalar med_p = `md'*`im'/(`total')*100
estadd scalar sobel_z = (`im' * `md') / sqrt((`im'^2 * `se_md'^2) + (`md'^2 * `se_im'^2))

local var Mexp_GMF_14after_std
local mv trust_govt
eststo m7: reghdfe y2     ///
`var'Xedu_nc     $cov##c.`var'    $agg##i.edu_nc  if `mv'!=.     ,a(country   i.birthyear) cluster(country)  keepsing
local total = _b[`var'Xedu_nc]
eststo m8: reghdfe `mv'     ///
`var'Xedu_nc  $cov##c.`var'   $agg##i.edu_nc  if `mv'!=.   ,a(country  i.birthyear) cluster(country)  keepsing
local im = _b[`var'Xedu_nc]  // stored for computing mediation effect
matrix SE = e(V)
local se_im = sqrt(SE[1,1])
eststo m9: reghdfe y2     ///
`var'Xedu_nc  `mv'  $cov##c.`var'   $agg##i.edu_nc    if `mv'!=.   ,a(country   i.birthyear) cluster(country) keepsing
local md = _b[`mv']
matrix SE = e(V)
local se_md = sqrt(SE[2,2])
estadd scalar med_p = `md'*`im'/(`total')*100
estadd scalar sobel_z = (`im' * `md') / sqrt((`im'^2 * `se_md'^2) + (`md'^2 * `se_im'^2))

estadd local cfe "√" , replace : m*
estadd local cxe "√" , replace : m*
estadd local mxi "√" , replace : m*

esttab m* using T3.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) /// 
keep (trust_govt Mexp_conVDem_14after_stdXedu_nc  Mexp_FOTPs_14after_stdXedu_nc Mexp_GMF_14after_stdXedu_nc) ///
varlabel(trust_govt "Mediator: trust in government" Mexp_conVDem_14after_stdXedu_nc "Media Freedom Exposure × No college") ///
mgroup("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}" , pattern(1 0 0 1 0 0 1 0 0)) ///
mtitles("{\i Trust in Science (IRT)}" "{\i Trust in government}" "{\i Trust in Science (IRT)}" "{\i Trust in Science (IRT)}" "{\i Trust in government}" "{\i Trust in Science (IRT)}" "{\i Trust in Science (IRT)}" "{\i Trust in government}" "{\i Trust in Science (IRT)}") ///
stats(med_p sobel_z cfe cxe mxi r2_a N, labels("Mediation effect as % of δ" "Sobel Z-value" "Country and birth year FE"  "Country-level controls × education" "Individual-level controls × Media freedom" "Adjusted R²"  "Observations") fmt( 2 2 0 0 0 3 0 ) )  ///
title({\b Fixed Effects Model: Mediation Effect of Political Trust}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Table 4. Mediation Effect of Media Trust
eststo clear

local var Mexp_conVDem_14after_std
local mv trust_journalists
eststo m1: reghdfe y2     ///
`var'Xedu_nc     $cov##c.`var'    $agg##i.edu_nc  if `mv'!=.     ,a(country   i.birthyear) cluster(country)  keepsing
local total = _b[`var'Xedu_nc]
eststo m2: reghdfe `mv'     ///
`var'Xedu_nc  $cov##c.`var'   $agg##i.edu_nc  if `mv'!=.   ,a(country  i.birthyear) cluster(country)  keepsing
local im = _b[`var'Xedu_nc]  // stored for computing mediation effect
matrix SE = e(V)
local se_im = sqrt(SE[1,1])
eststo m3: reghdfe y2     ///
`var'Xedu_nc  `mv'  $cov##c.`var'   $agg##i.edu_nc    if `mv'!=.   ,a(country   i.birthyear) cluster(country) keepsing
local md = _b[`mv']
matrix SE = e(V)
local se_md = sqrt(SE[2,2])
estadd scalar med_p = `md'*`im'/(`total')*100
estadd scalar sobel_z = (`im' * `md') / sqrt((`im'^2 * `se_md'^2) + (`md'^2 * `se_im'^2))

local var Mexp_FOTPs_14after_std
local mv trust_journalists
eststo m4: reghdfe y2     ///
`var'Xedu_nc     $cov##c.`var'    $agg##i.edu_nc  if `mv'!=.     ,a(country   i.birthyear) cluster(country)  keepsing
local total = _b[`var'Xedu_nc]
eststo m5: reghdfe `mv'     ///
`var'Xedu_nc  $cov##c.`var'   $agg##i.edu_nc  if `mv'!=.   ,a(country  i.birthyear) cluster(country)  keepsing
local im = _b[`var'Xedu_nc]  // stored for computing mediation effect
matrix SE = e(V)
local se_im = sqrt(SE[1,1])
eststo m6: reghdfe y2     ///
`var'Xedu_nc  `mv'  $cov##c.`var'   $agg##i.edu_nc    if `mv'!=.   ,a(country   i.birthyear) cluster(country) keepsing
local md = _b[`mv']
matrix SE = e(V)
local se_md = sqrt(SE[2,2])
estadd scalar med_p = `md'*`im'/(`total')*100
estadd scalar sobel_z = (`im' * `md') / sqrt((`im'^2 * `se_md'^2) + (`md'^2 * `se_im'^2))

local var Mexp_GMF_14after_std
local mv trust_journalists
eststo m7: reghdfe y2     ///
`var'Xedu_nc     $cov##c.`var'    $agg##i.edu_nc  if `mv'!=.     ,a(country   i.birthyear) cluster(country)  keepsing
local total = _b[`var'Xedu_nc]
eststo m8: reghdfe `mv'     ///
`var'Xedu_nc  $cov##c.`var'   $agg##i.edu_nc  if `mv'!=.   ,a(country  i.birthyear) cluster(country)  keepsing
local im = _b[`var'Xedu_nc]  // stored for computing mediation effect
matrix SE = e(V)
local se_im = sqrt(SE[1,1])
eststo m9: reghdfe y2     ///
`var'Xedu_nc  `mv'  $cov##c.`var'   $agg##i.edu_nc    if `mv'!=.   ,a(country   i.birthyear) cluster(country) keepsing
local md = _b[`mv']
matrix SE = e(V)
local se_md = sqrt(SE[2,2])
estadd scalar med_p = `md'*`im'/(`total')*100
estadd scalar sobel_z = (`im' * `md') / sqrt((`im'^2 * `se_md'^2) + (`md'^2 * `se_im'^2))

estadd local cfe "√" , replace : m*
estadd local cxe "√" , replace : m*
estadd local mxi "√" , replace : m*

esttab m* using T4.rtf , b(3) se(3) nonote nobaselevel  unstack noomitted  eqlabel(none) label  replace  star(+ 0.1 * 0.05 ** 0.01 *** 0.01) /// 
keep (trust_journalists Mexp_conVDem_14after_stdXedu_nc Mexp_FOTPs_14after_stdXedu_nc Mexp_GMF_14after_stdXedu_nc) ///
varlabel(trust_journalists "Mediator: trust in media" Mexp_conVDem_14after_stdXedu_nc "Media Freedom Exposure × No college") ///
mgroup("{\i Continuous Media Freedom Exposure (V-Dem)}" "{\i Press Freedom Exposure (FOTP)}" "{\i Media Freedom Exposure (GMF)}" , pattern(1 0 0 1 0 0 1 0 0)) ///
mtitles("{\i Trust in Science (IRT)}" "{\i Trust in media}" "{\i Trust in Science (IRT)}" "{\i Trust in Science (IRT)}" "{\i Trust in media}" "{\i Trust in Science (IRT)}" "{\i Trust in Science (IRT)}" "{\i Trust in media}" "{\i Trust in Science (IRT)}") ///
stats(med_p sobel_z cfe cxe mxi r2_a N, labels("Mediation effect as % of δ" "Sobel Z-value" "Country and birth year FE"  "Country-level controls × education" "Individual-level controls × Media freedom" "Adjusted R²"  "Observations") fmt( 2 2 0 0 0 3 0 ) )  ///
title({\b Fixed Effects Model: Mediation Effect of Media Trust}) ///
note("{\i \b Notes:} + {\i p} < 0.1, * {\i p} < 0.05, ** {\i p} < 0.01,*** {\i p} < 0.001 (two-tailed test)." )

**#Figure 2. Multilevel Model: Marginal Effect of Media Freedom Exposure on Trust in Science as Moderated by Science Literacy
local var Mexp_conVDem_14after_std
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.science_literacy || country:, mle
margins, dydx(`var') at(science_literacy=(-2.3(0.1)1.7))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("Continuous Media Freedom Exposure (V-Dem)") ytitle("Magrginal Effect of Media Freedom Exposure") xtitle("Moderator: Science Literacy") xsize(6) ysize(4) yscale(titlegap(-6))
graph save F2_mod1.gph, replace 

local var Mexp_binVDem_14after_std
mixed y2 female emp income loggdppc logpop log_nobel_pc log_qs500_pc e_pelifeex c.`var'##c.science_literacy || country:, mle
margins, dydx(`var') at(science_literacy=(-2.3(0.1)1.7))
marginsplot,recast(line) recastci(rarea) ciopts(fcolor(%30) lcolor(%1)) nolabels ///
title("Binary Media Freedom Exposure (V-Dem)") ytitle("Magrginal Effect of Media Freedom Exposure") xtitle("Moderator: Science Literacy") xsize(6) ysize(4) yscale(titlegap(-6))
graph save F2_mod2.gph, replace 

graph combine F2_mod1.gph F2_mod2.gph, ///
xsize(12) ysize(4) cols(2)
graph export "F2.png", replace
